Microscopic theory of non-adiabatic response in real and imaginary time 
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We present a general approach to describe slowly driven quantum systems both in real and imagi- 
nary time. We highlight many similarities, qualitative and quantitative, between real and imaginary 
time evolution. We discuss how the metric tensor and the Berry curvature can be extracted Irom 
both real and imaginary time simulations as a response of physical observables. For quenches ending 
at or near the quantum critical point, we show the utility of the scaling theory for detecting the 
location of the quantum critical point by comparing sweeps at different velocities. We briefly discuss 
the universal relaxation to equilibrium of systems after a quench. We finally review recent develop- 
ments of quantum Monte Carlo methods for studying imaginary-time evolution. We illustrate our 
findings with explicit calculations using the transverse field Ising model in one dimension. 
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I. INTRODUCTION 

Realizing efficient calculations of dynamical properties 
of interacting quantum systems remains one of the un- 
resolved challenges of modern physics. Even with some 
recent progress in simulating the dynamics in one dimen- 
sion using DMRG and related methods 1], as well as 
exact diagonalization (see, e.g., Ref. j2]), most physical 
systems remain currently out of reach. In a recent work 
[3] [5Q], we demonstrated that many difficulties can be 
overcome by going to imaginary time, where powerful 
quantum Monte Carlo (QMC) techniques can be used 
for a wide range of systems. This allows one to study 
some generic aspects of non-equilibrium dynamics and 
extract valuable qualitative and quantitative information 
pertaining also to real-time evolution of interacting sys- 
tems. We argued that the class of systems which can 
be analyzed in non-equilibrium setups therefore coincides 
with those which can be analyzed in equilibrium — those 
for which the QMC sign problems can be circumvented. 

Apart from numerical convenience, imaginary time dy- 
namics also has numerous experimental applications. If 
we are interested in the dynamics of a subset of degrees 
of freedom of a system, which couple to the rest of the 
system forming the environment, then the dynamics be- 
comes dissipative. There is no unique framework describ- 
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ing dissipative systems, but in many situations one can 
rely on Langevin dynamics, which is also equivalent to 
the model A dynamics [I] , which is in many cases equiva- 
lent to the imaginary time quantum dynamics with extra 
noise. If the bath has a temperature much lower than 
the driven system then the noise term becomes unim- 
portant. Then the model A dynamics describing the 
evolution of the real multi-component order parameter 
becomes equivalent to the imaginary time Schrodinger 
equation [4]: 



, dF 

dip. 



(1) 



where ipj is generally a multicomponent order parame- 
ter, the index j can be either discrete or continuous e.g. 
corresponding to the spatial coordinate, and F is the 
free energy of the system. This free energy can explic- 
itly depend on time if we are driving, e.g., external fields 
which explicitly enter F. Such situations were recently 
considered in Ref. [5]. Another wide range of applica- 
tions of imaginary-time quantum dynamics comes from 
applications to the Kardar-Parisi-Zhang (KPZ) equation 
and equivalent nonlinear Burgers equations which de- 
scribe the equilibrium behavior of polymers in random 
media, crystal growth, superconducting flux lines and 
many other systems (see, e.g., Refs. [HIE])- F° r instance, 
the differential equation describing the partition func- 
tion of a polymer in a disordered media takes the form 
of the imaginary time Schrodinger equation in a random 
potential, where the role of time is played by the coordi- 
nate along the polymer [7J. Using the replica trick, the 
KPZ equation maps to the imaginary-time Schrodinger 
equation describing bosons with short-range attractive 
interactions 8J. There are several other applications of 
Eq. ([I]) . In this work we will analyze the general proper- 
ties of the response of systems described by this equation 
together with the real time Schrodinger equation in situ- 
ations where the parameters of the system change slowly 
in time. 

The main purposes of this work are: (i) to further elab- 
orate our earlier findings of Ref. [3J , (ii) to give a quanti- 
tative comparison between real and imaginary time evo- 
lution for the specific case of the transverse field Ising 
model, (iii) to discuss how one can extract quantitative 
information about real time correlation functions from 
the low velocity asymptotics of the imaginary time re- 
sponse. We will also describe how one can extract real 
and imaginary components of the geometric tensor defin- 
ing the Riemannian metric and the Berry curvature asso- 
ciated with the ground state wave function from both real 
and imaginary time dynamics. We will discuss the appli- 
cation of the nonequilibrium scaling relations for quan- 
tum critical systems obtained in Refs. [3J HJ [TU] to 
accurately locate the quantum critical point and deter- 
mine the static and dynamic critical exponents from the 
collapse of physical observables. Our results can be use- 
ful for: understanding quantum annealing (in particular 
for finding the optimal path in the parameter space), ex- 



tracting long-time correlation functions and the dynami- 
cal exponent for disordered systems, evaluating both real 
and imaginary parts of the geometric tensor for inter- 
acting systems, including the Berry curvature and the 
fidelity susceptibilities. 

The paper is organized as follows. In Sec.[TT]we present 
the general theory of Kubo response of systems driven 
with constant velocity both in real and imaginary times. 
We identify the linear and quadratic susceptibilities of 
the response of arbitrary observables with respect to the 
velocity using Adiabatic Perturbation Theory, from here 
on abbreviated as APT. The linear susceptibilities are 
given by the components of the geometric tensor. In 
Sec. |III| we formulate the scaling theory for slowly driven 
gapless systems and systems driven through quantum- 
critical points and discuss its potential implications for 
experiments. We also relate the scaling of generic observ- 
ables to the scaling dimension of the relevant components 
of the geometric tensor. In Sec. |IV| we briefly review 
two complementary Monte-Carlo algorithms [3J [TT] for 
computing the quantities discussed above. In Sec. M we 
present the exact solution for imaginary-time quenches of 
the transverse-field Ising model in one dimension. From 
this solution we extract the scaling behaviors for several 
observables (e.g., excess heat, log fidelity, and nearest- 
neighbour spin-spin correlation functions). We compare 
these exact results with those obtained using APT and 
find a very good agreement between them. We also con- 



firm the general scaling relations presented in Sec. Ill In 



Sec. VI we compare the expectation values of different 
observables in imaginary versus real time quenches. In 
particular, for diagonal observables like the energy and 
the fidelity, we find a very good qualitative and quanti- 
tative agreement between the real and imaginary time 
evolution. For off-diagonal observables the low veloc- 
ity asymptotics in real and imaginary time are different 
(given by the imaginary and real components of the ge- 
ometric tensor) while at high velocities they only differ 
by a numerical factor of the order of one. In Sec. |VII| 
we illustrate how this scaling analysis, either in real or 
imaginary time, can be used to detect the position of 
the quantum-critical point without the need to vary the 
system size or the temperature. In Sec. |VIII| we discuss 
the universal behavior of the relaxation of observables 
(again in both real and imaginary times) to a prether- 
malized state following a quench near a quantum critical 
point. 



II. GENERALIZED KUBO RESPONSE OF 
DRIVEN SYSTEMS 

Originally Kubo response theory was developed to de- 
scribe transport coefficients for systems in weak electric 
fields. By now it refers to a general linear response theory 
to a static or time dependent external perturbation [12j . 
Let us point that in the Weyl gauge the scalar potential 
is zero and the electric field can be thought as the rate 
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of change of the vector potential: E = 1/cdtA. Thus, 
we can formally view the response to the electric field as 
the response to the rate of change of the vector potential. 
Here we extend this analogy and generalize the notion of 
Kubo response of a system to the rate of change of an 
arbitrary coupling A. 

Let us consider a system described by a Hamiltonian 
%[A(t)], where t is real or imaginary time and A is the 
vector of coupling constants in the parameter space. To 
distinguish between real and imaginary time we will re- 
serve the symbol r for the latter. In the real time case, 
the dynamics of the system is given by the time depen- 
dent Schrodinger equation: 

ia t m=nkt))m, (2) 

in imaginary time the dynamics is described by the cor- 
responding equation: 

d^(T) = -H(X(r))tKr). (3) 

We will assume that initially the system is prepared at 
t = or t = in the ground state. In imaginary time 
this assumption is not important since, for an evolution 
lasting sufficiently long time, the initial conditions be- 
come irrelevant. In real time our results can be readily 
generalized to arbitrary stationary initial conditions by 
performing statistical average of the expectation value of 
an observable over the adiabatically evolved initial den- 
sity matrix. We will also assume that the rate of change 
of the coupling is sufficiently slow, such that the system 
remains close to the instantaneous ground state (appro- 
priately rotated with the Berry phase) at all times during 
the evolution. We will also make the assumption that the 
ground state is not degenerate. In this case we can solve 
the Schrodinger equation using the APT. The details of 
the derivation of the first order corrections are presented 
elsewhere (see Refs. [HI [14] for real time and Ref. [3] for 
imaginary time). Here we extend the derivations to the 
second-order terms. 



A. Adiabatic perturbation theory in real time 

At the first order of APT in real time we find that the 
transition amplitude to the instantaneous state \n) 7^ |0) 
is given by |14) : 

«#>(*) = - / d*i(n|d tl |0)ex P H$ n o(ti,t)], (4) 
Jo 

where $ n o(ti,t) = $„(ti,t) — $o(^i, t) is the total phase 
difference accumulated between the ground and the ex- 



cited states in the time interval {t\ , i) , 

$„(*!,*) = / dt 2 (£ n (t 2 )-i(n\d t2 \n}), (5) 

where £ n (t) is the energy of the eigenstate \n) at time t. 
We emphasize that there is a difference in the limits of 
integration in Eq. ^ and Eqs. (12) and (13) in Ref. [14]. 
This difference is due to additional phase transformation 
in Eqs. (6) and (14) in Ref. Q3]. Going back to the 
original basis at the end of the calculation gives the result 
above. Note that the overall phase entering in Eq. ^ can 
be thought of as a purely dynamical phase coming from 
the gauge invariant energy: 

E n = £ n -v a A ( £\ (6) 

where v a = \ a and = i{n\d\ a \n) is the Berry con- 
nection associated with the n-th energy level. It is easy 
to see that E n is invariant under the gauge transfor- 
mation of the basis states by an arbitrary phase factor: 
|n) -> exp[i<f) n (X(t))]\n). 

Similarly in the second order of APT one finds 

4 2) (*)= Y, I dh f 1 dt 2 (n\d tl \m)(m\d t2 \0) 

x e -i[*„m(tl,t)+*mo(t 2 ,t)]_ (7) 

Let us point out that the first and second order terms 
in APT are not necessarily analytic functions of the adi- 
abatic parameter v = d t X. Thus, this is a non pertur- 
bative expansion. As an example of this, in Sec. [V] we 
will compare the results obtained within the first-order 
of APT with results from exact diagonalization and find 
a very good agreement, even when the observables have 
non-analytic dependence on velocity. In this sense APT 
does not have a formal expansion parameter, it only re- 
lies on the fact that the transition probabilities are small. 
For example, for a Landau-Zener problem the first order 
of APT gives the correct non-analytic dependence of the 
transition probability on the sweep rate, but leads to a 
small 7r 2 /9 — 1 deviation in the prcfactor [14] , 

For gapped systems, or for gapless systems in suffi- 
ciently high dimensions, the leading non-adiabatic cor- 
rections to various observables are analytic functions of 
the quench velocity |15j . Then it is possible to expand 
the expressions for the transition probabilities as a Tay- 
lor series in the velocity. This can be done noting that 
the integrand in Eq. Q is a product of a slow function 
(matrix clement) and fast function (phase factor) and in- 
tegrating by parts (see Ref. Q3] for details). Keeping 
terms up to the velocity squared we find that 



(n\d a \0) 1 d (n\dp\0) , . (n\d a \Q)(i{n\dp\n)-i{Q\d p \0)) 

p c v aV0 c W- nT" ~r -?T + ™ a V{, — — 
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where all matrix elements and energies are evaluated at 
time t. It is straightforward to check that if £ n =/= £ m 
then 



(n\d a \m) 



(n\d a H\m) 
f _ f 



(9) 



In Eq. ([8]) we neglected the additional fast oscillating 
terms which contain the initial excitations of the system. 
These terms can be suppressed either (i) if the protocol 
starts smoothly with zero rate or (ii) if the gap in the 
initial state is very large. The oscillating terms can be 
further suppressed because of various dephasing mecha- 
nisms. 

Similarly to the treatment above we can evaluate the 

(2) 

leading order contribution to ah. ■ Again, by neglecting 
the oscillating terms due to the initial excitations in the 
system, we find that for n ^ 0; 



,(2) 



(n\d a \m)(m\dp\0) 
(£ n — £o)(£ m — £q) 



(10) 



and finally the quadratic correction for the ground state 
amplitude reads: 



,(2) 



a n = - 



^1 
2 ^ I 

m/0 



,(D 



iv a Vp ^ 



(0\d a \m)(m\dp\0) 



Combining all the terms up to the second order in v we 
obtain the result: 



(n\d a \0) 



1 



d 



in 



V a Vp 



l^|0) 



(n\d a \0)(0\d ff \0)) 

- VaV " (£ n -£ r — VaV "^ (£ 



£n ~ £o d\ a £ n — £q 

v-^ (n\d a \m){m\dp\Q) 



- £o)(£m — £o) 



(12) 



where A nm (r) = £ n (r) - £ m {r) and 



a„{r) = o„(r)exp 



p dr'£ n (r') 



(14) 



This equation should be supplemented by the normaliza- 
tion condition: 



(15) 



at t = Tf, where t/ is the arbitrary final time of interest. 
From Eq. (14) it is clear that the coefficient a n (r) co- 
incides with the amplitude a n (r) for the evolved wave 
function to be in the instantaneous state \n) only at 
t = Tf. This boundary condition applies to the situ- 
ation where the dynamical process started in a distant 
past enough for the wave function to become insensitive 
to the actual initial state (which is always possible to sat- 
isfy in imaginary time). Note that, unlike the real time 
case, the integral equation ( 13 1 explicitly contains the 



unknown amplitude a n (r/), which has to be found from 
the asymptotic boundary condition. As in the real time 
case, if we deal with the eigenstates with a non-zero Berry 
connection, then one has to use the shifted (complex) en- 
ergies E n — £ n + v a (n\d a \n) = £ n — iv a A^ . Note that 
the fact that the energies are complex and the "moving 
Hamiltonian" is non-Hermitean is the consequence of the 
imaginary time evolution. 



From the integral equation ( 13 ) we find that in the 
leading order of APT: 



(i) 



(Tf) 



Tf 



dT(n\d T \0)e 



■P dr'A„ (r') 



(16) 



B. Adiabatic perturbation theory in imaginary 
time 



The APT analysis of the imaginary-time dynamics is 
very similar to that for the real time case. In Ref. [3] we 
derived the following exact integral equation for the am- 
plitudes of the wave function in the instantaneous basis: 

ot n {r) = a n (T f ) 

+ E P dT'{n\d T ,\m)c^(T')e-C f dT '' A ™( T "\ (13) 



For n^O and = 0, the first-order correction to the 
ground state amplitude vanishes. Similarly, in the second 
order of APT we find: 



,(2) 



(r/) = - £ JdT{n\d T \m)a^{T)e~^ dT ' 



A„ m (r') 



(17) 



As in the real-time case one can expand Eqs. ( 16 ) and 



(171 into a Taylor series in the quench rate: 



a 



(i) 



HjyO) 1 d (4gg|0) (n\d Q \0)((n\d p \n) ~ (0\d p \0)) (m s 

I 



where all energies and matrix elements are evaluated at t = Tf. Note that there is a sign difference compared 
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to Ref. [3J, due to a different definition of the velocity, 
v = —A, used in that work. Similarly, from Eqs. (13), 
([18} and (fl7| we find that for n ^ 0: 



(£„ — £o){£ m — £o) 



spin-spin correlation function, and so on (later on, when 
presenting the observables for the transverse field Ising 
model, we will introduce the following generalized forces: 
excess energy, the log-fidelity and the transverse magne- 
tization, those will be defined according to the definition 
in Eq. ( 22 1 in Sec. V B ). In the first two orders of the adi- 



Combining the two expressions above we find that up 
to v 2 terms and for n^O: 



abatic perturbation theory (both in real and imaginary 
time) we find: 



<*n(Tf] 



Md«|0) 



1 



() 



V a Vf3 



i\dp\0) 



My PS 



n\d a \m)(m\dp\0) 



<»ia |o)<oiaff|o) u y 

{£ n — Sq) 2 ~f (£ n — £$){£ m — £q 



E[K ( i 1) )*(^l^|0)+ai 1 )(0|a 7 H|n) 
E[(«l 2) )*("|9 7 H|0)+a( 2 )(0|9 7 H|n) 

n 

E (aPya£Hn\d y H\m), (23) 

We point that the real time expression for the transition where M (o) = _, \» H \0) is the ground state expecta- 
amphtude (12) can be formally obtained from the imagi- ^ yalue ^ & mQre gener&1 situation of finite initial 



(20) 



nary time expression above by the analytic continuation 
of the velocity to the complex plane v —> —iv. We ex- 
pect that this will be the case in all orders of expansion 
in the velocity. However, this continuation might not 
hold in general when there are additional non-analytic 
contributions like e.g. an exponential dependence of the 
transition amplitude on the velocity in the Landau-Zener 
sweep. There are no analogues of such exponential terms 
in imaginary time dynamics. We also point that in order 
to obtain the complex conjugate of the transition ampli- 
tude ( 12 ) from the imaginary time value ( 20 1 one needs to 



analytically continue velocity to positive imaginary axis: 
v — > iv. 

The correction to the amplitude for the n = state 
can be found by enforcing the normalization condition 
( 15 1 at t = Tf. 



(2) 



■-El 

2 ^ I 



,(D 



(21) 



C. Kubo response in the parameter space 

Having derived the expressions for the transition am- 
plitudes in real and imaginary time we can next compute 
the response functions. As in Refs. [31[TS], and without 
loss of generality, we will represent an observable as a 
generalized force: 



Mv 



(22) 



where 7 is some parameter in the Hamiltonian. The nor- 
malization factor in denominator highlights that in imag- 
inary time the wave function should be properly normal- 
ized. Clearly M 7 is simply the expectation value of the 
operator M. 1 = —d 1 H. For example if 7 is the exter- 
nal magnetic field, then the generalized force M 1 is the 
magnetization, if 7 is the volume then we have the pres- 
sure, if 7 is the spin-spin interaction then we have the 



temperature and real time dynamics, stands for 

the adiabatic expectation value of —dj'H, i.e., the ex- 
pectation value with respect to the density matrix adia- 
batically connected to the initial state. Combining this 
equation with Eqs. J8| and (10) for the real time case 
and with ( 18 1 and ( 19 1 for the imaginary case, we find: 



-in real time: 

M 1 « M 7 °) + F ja v a + [n 7Q;3 + I^ ap ]v a v p , (24) 



while in imaginary time we obtain: 



M 7 « M 7 °) - 2g~ fa v a + [ITL, 



ya/3 — n 7ct(3 JW Q W ( g. 



(25) 



To the linear order in the velocity, Eq. ( 24 ) was derived 



in Refs. [TBI EZJi an d in imaginary time it was derived 
in Ref. [3J. Note again the sign difference in the first 
term of Eq. ( [25] ) and the result in Ref. [3] due to dif- 
ferent sign conventions in the definition of v a . In this 
work v a = d T X a . In the above equations F ja and g ja 
are respectively the Berry curvature and the Riemannian 
metric tensor, which are related to the imaginary (anti- 
symmetric) and real (symmetric) parts of the geometric 
tensor [18] . Defining the geometric tensor as: 



we have: 



X«f) = (0fo|0) - (0|£;|0><0|fy|0), (26) 



F af 3 = i (Xap - Xp<x) = -23[Xa/9]> (27) 



9afS = 2 (Xafi + XP<*) = ^[Xafi) 



(28) 



The Berry curvature can be also expressed as a curl of 
the Berry connection [!§] : 

F a0 = d a A fj - d A a , A a = i(0\d a \0). (29) 

In the general case of a finite temperature, the expecta- 
tion values with respect to the ground state in the defi- 
nition of F and g should be substituted by the trace over 
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the density matrix representing the state adiabatically 
connected to the initial state. 

To the second order in the velocity the coefficients of 
the response functions are defined as follows: 



TT 1 



(0\d a \n)(n\d^n\m)(m\dp\0) 



(0|3 7 ^|0>]T 



(£„ — £o){£ m — So) 
(0\d a \n)(n\dp\Q) 



n^O 



(30) 



n 



y (0\d 7 H\n) d (n\d p \0) | ^ | 



y, (o|a 7 |n)(n|a a |o)(o|^|o) + r r 



<?ri — £o 



(31) 



It is straightforward to see that all the response coeffi- 
cients g ia , Frya , n^ Q/3 and 11^^ are gauge invariant, i.e. 
invariant under arbitrary basis transformations: 



MA)) 



There is a special class of observables in this regard: 
those that commute with the Hamiltonian in the final 
state. Examples of these observables include the Hamil- 
tonian itself (i.e., the energy of the system) as well as its 
various moments; the expectation value of any other con- 
served quantity; the probability to remain in the ground 
state, which is known as the fidelity, or in any particular 
eigenstate of the final Hamiltonian; the diagonal entropy 
of the system, and others. For these diagonal observables 
the linear response term in Eqs. (24) and (25) vanishes, 



as so does the response coefficient n^ Qj g. Therefore, the 
leading-order non-adiabatic response is given by the fol- 
lowing quadratic response function, which is the same 
both for real and imaginary time dynamics: 



M 7 °) 



(32) 



n^O 



(0\d a \n)(n\dp\0) 

fc ^~yz [\°irL)nn ~ {0 7 H)oo\ , 

(C n — CQ ) 



where (9 7 "H)„„ = {n\d 7 %\ 



(33) 



D. Relation to the zero-frequency limit of the 
response functions 

The components of the geometric tensor [20] and 
the second-order susceptibility II ' o can be expressed 
through the non-equal time correlation functions of the 
generalized forces. Using Eq. ([9]), the geometric tensor 
( 26 1 for a non-degenerate ground state can be rewritten 



n#0 



{0\d a H\n)(n\dpH\0) 



(£ n - £o) 2 



Next we recall the identity: 

i = _ lim r 

(£ n -£ ) 2 e^o+J 



-e£-»(e„-£o)S 



(34) 



(35) 



with which one can rewrite the geometric tensor ( 34 ) as 
an integral from the retarded correlation function: 



Xa/3 



dU(0\d a H(OdpH(0)\0) c e 



where the subindex c implies connected and 



(36) 



(37) 



is the Heisenberg representation of the generalized force 
d a T-L at the point of measurement. Note that here £ is an 
auxiliary variable, which is not related to the time evolu- 
tion during the dynamical process. One can also rewrite 
Eq. (36) in the imaginary time Heisenberg representation 
by formally rotating to Euclidean time: £ —¥ —it. 



One can also rewrite Eq. ( 36 ) through a derivative of 



the Fourier transform of the non-equal time correlation 
function: 



where 
G a p(uj) 



Xa/3 = -j9 w G Q/3 (w)| Lj=0 , (38) 



dte wt {Q\d a rl{t)d P rl{Q)\Q) c e- et . (39) 



From this expression we see that the metric tensor and 
the Berry curvature are given by the imaginary and real 
parts of the frequency derivative of the corresponding 
correlation function: 



(40) 



The components of the quadratic susceptibility II 7Q(9 
can be represented through the time-time correlation 
functions in a similar fashion. For example: 



n 1 

11 7 a^ 



/>oo />oo 

= / dti dt 2 e- £ ( tl+t2 ) t x t 2 
Jo Jo 

[{0\d a H(ti)a r Hd f ,'H(t2)\0) 

-<o|fl a ^(ii)a^(* a )|o)<o|a 7 ^|o)]. (4i) 



E. Equivalent derivation of the Kubo response via 
a generalized Galilean transformation 



The linear Kubo response given by Eqs. (|24| and (25) 



can be derived in a simple and intuitive way by going to a 
moving frame. The time dependent Schrodinger equation 



7 



^ can be rewritten in a comoving basis of the instanta- 
neous Hamiltonian as 



id t \ip) 



(H 



p 



(42) 



where the time derivative dt acts only on the coefficients 
of the expansion of the wave function in the comoving 
basis, while P a is a generalized momentum operator with 
respect to the parameter A Q . ft can be formally defined 
through the matrix elements in the instantaneous basis; 



(n\P a \m) — i(n\d a \m) = —% 



(n\d a H\m) 



(43) 



It is easy to see that P a is a Hermitian operator. This 
follows e.g. from differentiating the identity 



= S n 



(44) 



with respect to a. Alternatively one can note that the 
states \m(X)} can be obtained from some fixed basis cor- 
responding to e.g. Aq by some unitary transformation: 



|m(A)) = U mn \n(Xo). 



(45) 



For a non-degenerate spectrum this unitary operator 
corresponds to the adiabatic evolution of the Hamilto- 
nian. Then in the moving frame the Hamiltonian in the 
Schrodinger equation will clearly acquire an extra correc- 
tion 



iU^dtU = -ivJJ-^dJJ = ~v a P a . 



(46) 



Clearly the momentum operator P a = ill 1 d a U is the 
same as above in Eq. ( 43 ) . 



The RHS of Eq. ( |42[ ) extends the conventional Galilean 
transformation of the Hamiltonian in the moving frame. 
Indeed let us assume that we have a system of iV inter- 
acting particles in an external potential, which depends 
on time via: 

V{xt, ...x N ,t) = V(x! - X(t), ...x N - X(t)), 

where X(t) is a time dependent vector defining the mov- 
ing frame where the potential is stationary. This vector 
X(t) can also denote a center of mass coordinate of an 
interacting system. Then Eq. ( 42 1 is indeed the Galilean 
transformation where: 

v = X, and P = pj , 

are the usual velocity and momentum operator. 

Close to the adiabatic limit the additional term in the 
effective Hamiltonian in Eq. ([42]) can be treated as a per- 
turbation. A very similar equation can be written in 
imaginary time with the (non-Hermitean) Hamiltonian 
in the moving frame: H c ff = H — iv a P a . At the initial 
moment of time we turn on the velocity and at the mo- 
ment of measurement we effectively turn it off. Indeed, 



we may view the result of the instantaneous measurement 
as a result of a sudden quench of the velocity going back 
to zero, which is equivalent to going back to the original 
lab frame. Thus, for a simple protocol where the velocity 
suddenly turns on at time ti — and the measurement is 
done at time tf our perturbation in the effective Hamil- 
tonian looks like a pulse (see Fig.[T]). Now the results of 





' V = -v a P a 




f 

/ 
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/ 

/ 

* 







tf 



FIG. 1: (Color online) Schematic representation of the per- 
turbation in the effective Hamiltonian in th e moving frame 
for the case of constant velocity [see Eq. ( 42 1] . At the initial 
time ti = the velocity instantaneously changes from zero 
to a finite value and at the point of measurement t — tf the 
velocity effectively drops back to zero reflecting the transfor- 
mation back to the lab frame. The shape of the pulse can be 
smoothen at the initial time by turning on the perturbation 
slowly (dashed line). 

the APT, e.g., Eq. Q, leading to the Kubo formulas in 
real and imaginary time are easily understood as those 
of the ordinary perturbation theory with respect to the 
velocity dependent term in the effective Hamiltonian; see 
Eq. (42). Indeed, if the time of the pulse is short, then 



we can apply an ordinary time dependent perturbation 
theory where the transition amplitude is proportional to 
the matrix element of the perturbation integrated over 
time, i.e.: 



v a P a dt ~ 6 \ a P a . 



(47) 



Thus, in this limit of a short pulse we simply reproduce 
the ordinary perturbation theory where the transition 
amplitude is proportional to the change of the coupling 
constant. Conversely, in the long time limit, the response 
is very different: to leading order of perturbation theory 
the system is excited at the initial time, then it freely 
evolves in a moving basis, after which it is excited again 
by the quench at t — tf. Thus, the transition ampli- 
tude is the sum of two terms, corresponding to the initial 
and final quenches. The contribution due to the initial 
quench comes in as a rapidly changing Rabi phase. As we 
discussed above, this phase averages to zero due to any 
dephasing mechanism, or it can be simply suppressed by 
a slow turning-on of the velocity. If this is the case, then 
we can reproduce the transition amplitudes (|8|) and ( 18 ) 



by using an ordinary static perturbation theory with re- 
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spect to the perturbation V = v a P a for real time proto- 
cols and V = iv a P a for imaginary time protocols. 



F. Extracting the Berry curvature from the 
imaginary time dynamics. Analytic continuation of 
the imaginary time dynamical response to real time. 



Equations (24 1 and (25) show that one can extract 



the Berry curvature as a linear response to the quench 
velocity in real time and the metric tensor from the 
imaginary time response. Real time dynamics of course 
corresponds to physical processes and thus can be di- 
rectly realized in experiments. However, with the ex- 
ception of one dimensional systems, there are very lim- 
ited numerical tools which would allow for real time sim- 
ulations. Imaginary-time dynamics has the advantage 
that it is amenable to powerful QMC as we first showed 
with the non-equilibrium QMC (NEQMC) approach in 
Ref. [3] and later with a different quasi-adiabatic QMC 
(QAQMC) scheme in Ref. [TTJ. We will briefly review 
these approaches in Sec. [TV] 

Because of the opportunities offered by the NEQMC 
and QAQMC methods, it would be very practical to have 
a way to extract the Berry curvature using imaginary- 
time dynamics as well. This indeed becomes possible 
if we evolve left and right wave functions in imaginary 
time with different velocities, i>„ and , and evaluate 
the following expectation value: 



i [(Mr)\dpH\^ R ( T )) - (^(r)|9^|^(r))] , 

where iPl(t) and iPr(t) are the solutions of the imaginary 
time Schrodinger equation for the two protocols charac- 
terized by different velocities. It is important that both 
tPl(t) and iPr(t) are evaluated at the same point in the 
parameter space. Then it is easy to see that in the linear 
order in the velocities and vr we have: 

« ^a(i£ - «£)• (49) 

In particular, if — » 0, i.e. if the left wave functions is 
the instantaneous ground state, then the equation above 
is equivalent to the real time linear response (24). Very 



similar results apply to real-time simulations, i.e., if we 
compute Mp in a real-time protocol we will get the lin- 
ear response proportional to the real part of the metric 
tensor. 



The result ( 49 ) opens the possibility of extracting the 



Berry curvature from numerical simulations in imaginary 
time. Note that the Berry curvature is non-zero only if 
the time reversal symmetry is broken, and, hence, the 
wave function is complex. In these situations QMC sim- 
ulations usually suffer from the sign problem. However, it 
is important to mention that the time reversal symmetry 
can be broken by the parameter /?, which does not enter 
the time evolution and only appears in the definition of 
the observable we measure. In these situations, the wave 
function always remains real and the sign problem may 
be avoidable. 



In imaginary time dynamics left and right states cor- 
responding to opposite velocities naturally occur for the 
asymmetric expectation values. For a closely related 
QAQMC algorithm this issue was discussed in detail in 
Ref. [TT] . Let us define the protocol where the coupling 
A changes in imaginary time in the symmetric fashion in 
the interval: [0,2T]: A(r) = A(2T-r). Then denoting 
the imaginary time evolution operator 



U(t 1 ,t 2 ) = T T exp 



U(T)di 



(50) 



we can write a generally asymmetric value of arbitrary 
observable M 7 (r, 2T — t) as 



M 7 (t,2T-t) = 



(^ \U(2T,t)M 7 U(t,0)\tI; ) 



(51) 



(Vo|tf(2T,0)|V>o) 

where ipo is the initial state and r G [0, 2T]. Such ex- 
pectation values are very straightforward to evaluate for 
any t in e.g. Monte-Carlo simulations (see also discus- 
sion below following Eq. (76)). It is easy to see that in 
the middle point of the evolution r = T this expectation 
value becomes equivalent to the result of the imaginary 
time dynamics discussed in the previous sections. Away 
from the symmetric point t =/= T the left and right states 
effectively evolve with the opposite velocities so the ex- 
pectation value above effectively becomes 

^ (M-v)\M 7 \Mv)) 



M 7 (r, 2T-r) 



(ip L (-v)\ip R (v)) 



(52) 



As it is discussed in Ref. [TT] the interval around r = T 
where the crossover between symmetric and asymmetric 
asymptotics happens vanishes as v — ¥ 0. 

From the analytic properties of the wave-function dis- 
cussed below Eq. ( 20 ) it is clear that the real time expec- 



tation value of the observable M 1 can be obtained by the 
analytic continuation of Eq. ( 52 ) to the imaginary veloc- 
ity v —¥ —iv. Indeed as we noted earlier to get the real 
time result for the wave function and its complex conju- 
gate one needs to analytically continue the velocity to the 
imaginary axis in the opposite ways v — > ±iv. However, 
because in Eq. (52 1 ^ and ipR are evaluated at opposite 
velocities both have to be analytically continued in the 
same way v — > —iv. As a result the analytic continuation 
works for the observable M 7 . It is interesting that for- 
mally this analytic continuation is valid perturbatively to 
all orders in v. In particular, this implies that the leading 
asymptotic of off-diagonal observables in the asymmetric 
points in (51 1 will be given by the Berry curvature mul- 



tiplied by the velocity and for the diagonal observables 
the leading asymptotic will be quadratic in velocity but 
will have a negative sign (i.e. opposite in sign to the real 
time asymptotic). 



G. Extension to non-linear quench protocols 

The linear response analysis carried out above was 
based on the assumption that near the point of mea- 
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surement the quench velocity v a is non-zero. This is the 
case for generic protocols. However, it is also possible 
to especially design protocols where one approaches the 
point of measurement t — tf with some other power law 
characterized by the exponent r > (see also Ref. 21): 



A(f)« Xf+Vr ^ t)r . 



Here the parameter: 



-l) r d r X{t)/dt r \ t=tf 



(53) 



(54) 



is the adiabaticity parameter, which plays the role of the 
quench amplitude for sudden quenches (r = 0), quench 
velocity for linear quenches (r = 1), quench acceleration 
for quadratic quenches (r = 2), and so on. We should 
note that, with the definition ( |53[ ) for a linear quench, 

v r = —A is actually a negative velocity. However, the 
convention above is more natural for non-linear protocols, 
in particular for noninteger values of r. In Ref. [5] it was 
also analyzed the case when r is negative and were found 
different behaviors. In this work we will not be concerned 
with such situations. 

To evaluate the leading non-adiabatic response for 
quenches according to (53) with arbitrary r we need to 
perform the asymptotic analysis of Eqs. Q and ( 16 1 us- 
ing Eq. (53) for the dependence of A(i) and taking the 
asymptotic limit \v r \ — > 0. In this limit the matrix el- 
ements and the energy spectrum entering Eqs. Q and 
(16) are approximately constant. Then, as an example, 



in imaginary time Eq. ( 16 ) reduces to: 



« = -(n\d a \0) / dX a e- A ^f- T ^\ 



(55) 



where the matrix element (n|9 a |0) and the energy dif- 
ference A n o are computed at the point of measurement: 
A = Xf. The integral is taken along the actual path A(r) 
and is readily evaluated in the long-time limit: 



(n\d Q \0) 



(56) 



For the linear quench this expression reduces to the lead- 
ing term in Eq (18), keeping in mind that for r = 1 we 



have v r = —v; see Eq. (53). Similarly, in real time we 
find 



( - i) r Vr,a 



(n|fl a |0) 



(57) 



Using these excitation amplitudes we can easily find lin- 
ear response expressions for the generalized forces ex- 
tending Eqs. p4j) and ((25}. 



In real time: 
while in imaginary time: 



i) r X r cX l ) 



v a + 0(v 2 ), (58) 



0(v 2 



(59) 



Here we have defined: 

(0\d y H\n)(n\d a H\0) 



X, 



;r = E 



(£ n - £ ) r+1 



t 



dt-e- et (0\d.yH(t)d a H(0)\0). (60) 



As for the linear quenches the difference between real 
and imaginary time transition amplitudes is contained 
in the phase factors. One can always eliminate these 
factors artificially by considering separately left and right 
states, e.g., in the imaginary time evolution and forming 
the appropriate linear combinations similarly to Sec. |II F| 
For diagonal observables such as the energy or generic 
observables which are measured not instantaneously after 
the quench at t = tf but after allowing the system to 
relax to the diagonal ensemble, the phases in Eqs. (58) 
and (59) do not matter. Thus, to leading order in v the 



responses in real and imaginary time coincide and are 
given by: 



M 7 = (V|A^ 7 |V) 



n 



(61) 



where: 



n 



(0\d a n\n)(n\dpH\0) 



£ ) 2r+2 



(n\M 7 \n) - (0|X 7 |0)] . 

(62) 



This expression opens a way of measuring the symmetric 
part of the geometric tensor in real-time experiments. 
For example, one could measure the excess heat: M. 1 — 

H for a protocol with r = 1/2. In this case n^,'^ 2 reduces 
to the metric tensor g a p. 

Furthermore, we note that if the power of the quench 
is r = 4, the equations (58) and (59) become identical. 



This suggests the interesting fact that by probing a sys- 
tem through a quartic quench (~ i 4 ), we get the same 
response in real and imaginary time dynamics. 



III. UNIVERSAL SCALING NEAR QUANTUM 
CRITICAL POINTS 

The linear response theory presented above allows one 
to associate deviations from adiabaticity of various ob- 
servables through different susceptibilities. These suscep- 
tibilities are in turn expressed through integrals of non- 
equal time correlation functions and, in particular, are 
very sensitive to their long time asymptotics. However, in 
gapless regimes, specifically, near continuous phase tran- 
sitions, these susceptibilities may diverge and the linear 
response theory then breaks down. In these situations 
one can extend standard scaling theory of continuous 
phase transitions to non-equilibrium setups. 

The scaling hypothesis originally introduced by 
Pokrovsky-Patashinsky and Kadanoff (see, e.g., Refs. [22j 
I23j for references) is based on the conjecture that univer- 
sal physics near continuous phase transitions depends on 
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the microscopic parameters only through the correlation 
length. For quantum phase transitions, which are rele- 
vant for our work here, the correlation length universally 
diverges with the tuning parameter |23j as: 



&(A) 



1 

ia-a c k 



(63) 



where v is the correlation length exponent. For a multi- 
dimensional parameter space one can have different ex- 
ponents v along different directions. In this case there is 
more than one correlation length. Then, in finite size sys- 
tems the scaling hypothesis states that the expectation 
value of any observable can be written in the following 
scaling form: 



M 7 = const + L-^f^L/txiX)), 



(64) 



where the first constant term represents a non-universal 
non-critical contribution to M 7 , while /i 7 is a universal 
number which defines the scaling dimension of the oper- 
ator A4 7 , L is the system size and / 7 (x) is a universal 
function. Various two-point correlation functions can be 
represented in a similar form |23] where instead of the 
system size we use the separation between points. By 
definition the scaling dimension of the tuning parameter 
is 1/v. 

In general, the scaling dimensions of the coupling A a 
and the corresponding generalized force M a are inde- 
pendent. But there is an important exception, namely, 
when the coupling X a is relevant, i.e, this coupling drives 
the system to or away from the critical point. Then 
the scaling dimension of the product (A — \ c )a-M a must 
be equal to z, the scaling dimension of the energy, and 
thus we must have A Q + 1/v = z. If we are dealing 
with marginal perturbations which keep the system in 
the gapless regime, e.g. which lead to renormalization 
of the velocity, then we generally have to set v — > oo 
(corresponding to scaling dimension of A Q equal to zero) 
leading to A Q — z. In this work we will assume that 
the driving term in the Hamiltonian is either relevant or 
marginal but that the observable M. 1 has an arbitrary 
scaling dimension /i 7 . The generalization to situations 
where the driving term is irrelevant is straightforward. 

As it was argued in Refs. [1 G2 MTT\ I3T] [M] the scal- 
ing ansatz ( 64 1 can be extended to non-equilibrium sit- 



uations if we add the quench velocity as another scaling 
variable. In particular, since v a = d\ a /dt we can expect 
that the scaling dimension of the velocity is dim[i; Q ] = 
dim[A Q ] — dim[t] = v + z. This scaling dimension im- 
plies that there should be another length scale associated 
with the velocity v: 



(65) 



For generic power law protocols characterized by the ex- 
ponent r [see Eq. (53 1] the above result immediately gen- 
eralizes to: 

6u ~ V"^+T. 



This length scale is indeed the Kibble-Zurek length, 
which was first introduced by Zurek |25j for classical 
phase transitions in relation to the Kibble-Zurek mecha- 
nism and later reintroduced in Refs. [37] HH] for 
quantum phase transitions. Physically describes the 
length scale beyond which the system effectively freezes 
and cannot follow the instantaneous ground state. For 
multicritical points with more than one direction corre- 
sponding to different critical exponents, £ v can depend 
on the direction of the quench [29] . Now one can use this 
additional length scale to extend the ansatz for both real 
and imaginary time evolution (64) to: 



M 7 (A, v) = const + L-^/ 7 (£/6(A), ( 66 ) 

Here A and v describe the coupling and its rate of change 
at the point where we perform the measurement of our 
observable M 7 . 

The asymptotics of the scaling function / 7 can be of- 
ten determined from qualitative considerations. Thus, if 
£u ^ £x, then the system effectively behaves adiabat- 
ically and we should recover the scaling behavior per- 
taining to the static equilibrium. In the opposite limit, 
depending on the ratio of and L, we should recover 
a similar crossover between the linear response discussed 
earlier and non-equilibrium universal scaling. For £ v 3> L 
we have fJL/£ v ) ~ (L/£ v ) z+1 /» cx v if there i 



is a non- 



vanishing relevant component of the geometric tensor or 
/ 7 (L/^) ~ (L/£ v ) 2z+2 < v cc v 2 in cases where we ex- 
pect quadratic scaling with the velocity, e.g., if M 7 is 
a diagonal observable. This asymptotic predicts a non- 
trivial scaling of the observables with the system size for 
very slow protocols. In the opposite limit L » we 
expect the extensive observables to scale linearly with 
the system size, i.e., f^(L/^ v ) ~ ) M " T+1 - While 

for intensive observables M 7 should saturate to a con- 
stant value independent on the system size, such that 

These simple considerations well reproduce the scaling 
behaviors derived earlier. For example, we expect that 
the density of defects is intensive and that the number of 
defects has scaling dimension zero (more accurately this 
statement applies to the log fidelity [21]) and thus: 



L d 



(67) 



is the well known Kibble-Zurek scaling form [35]. Simi- 
larly we can recover the scaling of the excess energy den- 
sity for quenches ending near a quantum-critical point: 
Q ~ v( d+z ^/( zu+1 \ which follows from noting that the 
scaling dimension of the Hamiltonian is z (or, equiv- 
alently, that the scaling dimension of the Hamiltonian 
density is d + z). These scaling considerations equally 
apply to imaginary-time (dissipative) and real-time dy- 
namics. The only difference is that in the low-velocity 
limit £.„ <C L the response is given by different suscep- 
tibilities, which, however, have the same scaling proper- 
ties. Note that while in this paper we focus on quantum- 



11 



critical points, the general considerations equally apply 
to thermal transitions, as it was emphasized in Ref. [5]. 



IV. TIME-EVOLVING QUANTUM MONTE 
CARLO ALGORITHMS 



One of the primary reasons for considering imaginary- 
time dynamics is that it is amenable to numerical simu- 
lation with modified QMC methods. This way, one can 
go beyond one dimension (where DMRG is applicable in 
real time) for a rather broad class of systems for which 
sign problems can be avoided. This class coincides with 
that for which equilibrium QMC methods can be applied. 

Standard QMC algorithms can be classified into finite- 
temperature methods, where the goal is to compute a 
quantum-mechanical thermal average of the form: 

(A) = ~Tc{e- pH }, Z = Te{e- pH }, (68) 

and ground-state projector methods, where some oper- 
ator P((3) is applied to a "trial state" l^o)) such that 
\~9ff) — P(/3)|* ) approaches the ground state (up to an 
irrelevant normalization factor) when j3 — > oo and an 
expectation value (which is what normally is computed, 
although one can also stochastically generate the wave 
function) 



(A) = ±(*p\A\*p) 



Z=(*p\* l ; 



(69) 



approaches the true ground state expectation value, 
(A) —> (0\A\0). For the projector, one can use P(/3) = 
e~' 3 ^ with large fi or a high power of the Hamiltonian, 
P(m) = H m . If in the latter case one uses m cx N/3 with 
/3 of the former approach, the same rate of convergence 
applies for a given system volume N (which follows, e.g., 
from a series expansion of the exponential, which for large 
(3 is dominated by powers of the order n — /3\Eq\, where 
E is the ground state energy). 

The differences between T > and T = projector 
methods can be thought of in terms of different boundary 
conditions in imaginary time: The trace taken at T > in 
( 68 1 corresponds to periodic boundaries while the projec- 



tor methods correspond to opening up these boundaries 
and replacing them with the ones corresponding to the 
trial state. Completely open boundary conditions corre- 
spond to the trial state being the equal superposition of 
all states in the basis used. 

The time-evolving QMC methods we have developed 
are essentially modified projector algorithms. In the orig- 
inal NEQMC approach the exponential operator e _/3w 
for a fixed Hamiltonian is replaced by the Schrodinger 
evolution operator in imaginary time, 



U(t) = T T exp 



dT'H[\(T')} 



(70) 



where T T indicates time ordering. As in equilibrium 
QMC schemes, there are several ways to deal with the 



exponential. In the context of spins and bosons, the most 
frequently used methods are based on (i) the Suzuki- 
Trotter-decomposition, which leads to world-line meth- 
ods |30j . (ii) the continuous-time version of world-lines 
(e.g., the worm algorithm |31j). and (iii) the Taylor ex- 
pansion leading to the SSE method [32l [33] (see Ref. [34] 
for a recent review of these approaches) . The latter two 
methods are not affected by any approximations (beyond 
statistical sampling errors) , while (i) has a discretization 
error. 

In the NEQMC algorithm a series expansion is em- 
ployed for (70 1 in the non-equilibrium case, while in the 



more recently introduced QAQMC method the power 
T~L m of the Hamiltonian used in standard projector meth- 
ods is replaced by a product of evolving Hamiltonians. It 
was shown in Ref. |llj that the product evolution repro- 
duces imaginary-time Schrodinger dynamics up to the 
leading corrections in v to the adiabatic evolution. In 
practice, this kind of method is easier to implement than 
the NEQMC scheme, and, moreover, one can obtain re- 
sults for all times between the initial and final Hamil- 
tonian in a single run. We here briefly review the two 
methods. 



A. Non-equilibrium Schrodinger approach 

In the NEQMC scheme first proposed in [3] the expo- 
nential in ( 70 1 is expanded in a power-series and applied 



to an initial state ^(O)): 



oo „ T 

i*w> = E / dr « 



dr n -i ■ ■ ■ dr 2 x 



[ 2 dr 1 [-n(r n )}---[-rl(n)}\m)- (71) 

J To 



'T 

The Hamiltonian is a sum of terms, 



(72) 



where the index i can refer to lattice sites, links, etc., 
and N op is the total number of these operators. A minus 
sign has been included for convenience. The operator 
product in ( 71 ) is then written as a sum over all strings 



of the operators Hi. Truncating at some maximum power 
n = M (adapted to cause no detectable error — see [S3] 
for a discussion of this issue in the SSE method) and 
introducing a trivial unit operator H = 1, one obtains: 



|*(r)> - 



(M-n)\ 



E ( T - To )M-n 



dr„ 



dr-2 x 



f * dr 1 H im (T m )---H i2 (T 2 )H il (T 1 )\y(0)), (73) 

J T 

where i p € {0, ...,M}, J2h denotes the sum over all 
possible sequences of M operators Hi and n is the number 
of indices i p ^ 0. 
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At this stage a basis {\a)} should be chosen. For spin 
systems, this would normally be the standard basis of 
the z spin components: |a) = \S*, . ■ . , S^). If the initial 
state (^(O)) has some simple known expansion in this 
basis, |^(0)) = J2a c a\ a )> tms can be used in (73). The 



scheme is particularly simple when using the equal su- 
perposition, e.g., | *(())) = Ili(ti + U) f° r an 5 = 1/2 
system, but other states can be used as well. For models 
with spin-isotropic interactions, such as the Heisenberg 
model, it is easy to use amplitude-product states in the 
singlet sector [35]. One can also start with the ground 
state of some Hamiltonian H(\q), by adding to (71) a 
projection with that fixed Hamiltonian before the quench 
with the time dependent H[X(t)] is applied. 

For practical implementations of the QMC scheme, the 
terms Hi of H should have the property that Hi\a) — 
hi(a)\a'), where \a') is a basis state. In the standard 
spin basis, this implies that H t is either a diagonal or 
off-diagonal operator (i.e., i denotes not only a lattice 
unit but also refers to either a diagonal or off-diagonal 
part). A string of operators and their associated time 
values, along with a state \a) then constitute a configu- 
ration, and the QMC simulation amounts to importance- 
sampling of these configurations, which strongly resemble 
those of a path integral. 

To guarantee the absence of a sign problem we need to 
place certain conditions (which are not always possible to 
satisfy) on the matrix elements hi(a), i.e., the product 



of all matrix elements corresponding to a term in ( 73 ) 



has to be positive. While this is a limitation of the QMC 
approach in general, the class of accessible models is still 
large and includes highly non-trivial and important sys- 
tems (see Ref. |36j for a recent review of quantum spin 
models without sign problems). 

Expectation values (t)\A\^ (t)} / (*(t)|*(t)) are 
computed by sampling the normalization (^>(t)\^(t)) 
written with ( 73 1 . The procedures are very similar to 



those used in the standard SSE and ground-state pro- 
jector methods [3H [35J [37] . The main difference is that 
each operator is associated with a time value. The sim- 
plest way to deal with the times is to sample them com- 
pletely independently of the operator and state updates, 
i.e., the latter are performed with fixed time values, ac- 
cording to one of the standard schemes [3U [35] [37] , and 
the times are updated without changes in the operators 
and states. A segment of m times, t,, . . . , T i+rn ^i, can be 
simultaneously updated by generating m numbers within 
the range (Tj_i, Tj_|_ TO ), then order these times according 
to a standard scheme scaling as log(m) [38]. The ordered 
set is then inserted in place of the old set of times, with 
a Metropolis acceptance probability obtained from (73 1 . 



The number m can be adjusted to give an acceptance 
probability close to 1/2. 

As discussed in Sec. |II F[ in addition to conventional 
expectation values, it is also useful to study asymmetric 
expectation values defined with two different time evo- 
lution operators U and V, e.g., corresponding to differ- 
ent velocities: (ipo\V*AU\ip )/(ipo\V*U\ij ). This can be 



done with a simple generalization of the above NEQMC 
algorithm. 



B. Quasi-adiabatic approach 

One may ask whether the role of the time integrals in 
Eq. ( 71 1 is crucial. Clearly, they are needed in order to 



obtain a mathematically exact expansion of the time evo- 



lution operator (70), but one can, in fact, also formulate 
a scheme similar to time evolution without these inte- 
grals, by acting on the initial ground state l^o) 01 %[Ao] 
with a product of M evolving Hamiltonians: 

Pali = [-tt(A M )]....[-ft(A 2 )][-ft(Ai)] > (74) 
where we consider the parameter changing according to 



A 4 = A + tA) 



(75) 



and A\ = [Xt+i — A$]/M is the single-step change in the 
tuning parameter. One can also consider a non-linear 
grid of points, but here we focus on the linear evolution. 
It is clear that |V>m) = Pm, i|*o) approaches the ground 
state I^m) of the final Hamiltonian W[Am] m the limit 
of large M (up to an irrelevant normalization). 

In Ref. [TT] it was also demonstrated, using APT, 
that NA\, where N is the system volume, plays the 
role of a velocity v, and that \i/jm) captures the lead- 
ing non-adiabatic corrections in v to the imaginary-time 
Schrodinger evolution. This is sufficient for recovering all 
the dynamical susceptibilities that we discussed in this 
work. 

Moreover, one can also consider generalized (asymmet- 
ric) expectation values of the form: 



(tt(A )|Pi,AfP A f,t+iAP M |tt(A )) 
(*(A )|Pi,m-Pm,i|*(A )) 



(76) 



where only the special case t = M corresponds to a 
true quantum mechanical expectation value but also the 
generic t ^ M quantities contain useful dynamic infor- 
mation and obey dynamic finite-size scaling. A signifi- 
cant advantage of QAQMC over the NEQMC approach 
is then that one can compute (A) t for all t simultaneously 
in a single simulation for operators A that are diagonal 
in the basis used. Such a simulation amounts to gen- 
erating terms (paths) contributing to the normalization 
( , J(Ao)|-Fi i m-Pm,i| , I , (Ao)) and successively measuring di- 
agonal observables after propagation of the state with t 
operators, for t on a suitable grid. 

Figure [2] shows examples of results obtained with the 
QAQMC method in simulations of the ID transverse-field 
Ising model, which we introduce in the next section. The 
quantity shown is the magnetization fluctuation, 



X = N((ml)-(\rn z \) 2 ) 



(77) 



which exhibits a peak close to the known quantum- 
critical point at g = J. The peak grows both as a func- 
tion of the size L and m, and one can subject the data 
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FIG. 2: (Color online) Magnetization fluctuation, Eq. (77 1, of 
the ID transverse-field Ising model (see Eq. (78)) in QAQMC 
calculations with different length M of the operator string for 
an L = 32 system (left) and for different system sizes at fixed 
M (right). The Ising and field terms are J = A, g — 1 — A, 
so that the quantum-critical point is at A = 1/2. The whole 
< A < 1 curve was obtained in a single simulation. 



to various forms of finite-size and finite-velocity scaling, 
examples of which are discussed in Ref. [TT] . 

As an alternative to computing expectation value 
based on evolving the same state from the left and the 
right in ( 76 ) , one can also carry out QAQMC simulations 



as a one-way evolution. In the simplest case, the left state 
(i/>l| is the ground state of A and the right state {if>n\ 
is the ground state of Am • The single sequence ( 74 ) be- 



tween these state will then smoothly connect them, and, 
again, this evolution captures the leading non-adiabatic 
corrections in v to the standard Schrodinger dynamics, 
with v oc LA\. As we discussed earlier (see Sec. I1F) 



asymmetric expectation values including one way evolu- 
tion can be used for computing the Berry curvature in 
the system. 



TRANSVERSE FIELD ISING MODEL IN 
ONE DIMENSION 



The ID transverse-field Ising model is defined by the 
Hamiltonian: 



H 



(78) 



where o~ x and o~ z are Pauli matrices and (ij) are nearest 
neighbours sites. The parameters g and J control the 
nature of the quantum state: for gj J > 1 the system 
is a quantum paramagnet, while for g/J < 1 it is in a 
magnetically ordered phase. The point (g/J) c = 1 of this 
spin chain is a quantum critical point (QCP) separating 
the two different phases. 



A. Imaginary time quench of the transverse field 
Ising model: exact solution 

We investigate the dynamical response of the system 
in the vicinity of the QCP by changing in time either 
g or J. Quenches in real time of the Ising model have 
been considered previously [14] . while less is known in 
imaginary time. In Ref. [3] we have fixed g = 1 and 
considered a particular quench protocol for J: J(t) = 1 + 
vt, starting in the ground state at tq = — 1/v and ending 
at the QCP at r = 0. We can implement a similar process 
fixing J = 1 and changing g as g(r) — 1 — vt. If instead 
we change g as g(r) = 1 + vt, then we approach the 
QCP from the initially ordered (ferromagnetic) phase. 
All these protocols give a very similar scaling behavior of 
the observables. Thus, for extracting analytical results 
we will focus on the particular protocol: 



9 (t) = 1 + X(t), 



A = -i 



(79) 



In the next section we will illustrate our results with nu- 
merical simulations, in which we also consider different 
protocols. 



We now investigate the results of the previous sections 
by analyzing quenches in real or imaginary time of the 
one-dimensional (ID) transverse field Ising model. This 
model maps onto free fermions and, thus, it is easily solv- 
able. It was used to rigorously demonstrate the univer- 
sal scaling relations both for equilibrium phase transi- 
tions [23] and various aspects of quantum dynamics, in- 
cluding the Kibble-Zurek scaling [55] . Furthermore, us- 
ing a closely related model it was recently demonstrated 
that the universality of slow quantum dynamics does not 
rely on integrability [TU]. We also point out that this 
model was extensively used to study the dynamics fol- 
lowing sudden quenches (see e.g. Refs. [40H43] ). Here we 
will use this model again for the purpose of a detailed 
comparisons between real and imaginary time dynamics, 
to establish that the universal aspects are identical. We 
will also demonstrate how one can use this universality 
to accurately extract the equilibrium transition point and 
the critical exponents using non-equilibrium protocols. 



Spectrum of the Ising chain 



It is well known that the Hamiltonian (78) can be 



mapped to that one of non-interacting fermions using the 
Jordan- Wigner transformation [23 , 39 . Because of trans- 
lational invariance, the relevant excited states are only 
those which contain pairs of quasi-particles with oppo- 
site momenta. As a result, in this reduced Hilbert space 
the Hamiltonian of the system splits into a direct sum 
of Hamiltonians describing two-level systems with the 
states: f)fc an d | i)k corresponding to empty and filled 
fermionic levels with momenta (k, —k), respectively: 



fc>0 



where: 



H k = -2[g - cos(fc)]o- z + 2sin(k)a x . 



(80) 



(81) 
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Each of these Hamiltonians has the following eigenvec- 
tors: 



2. Linearized spectrum: exact solution 



sin(f? fc /2) 
-cos(6> fe /2) 



cos(6» fe /2) 
sin(0 fe /2) 



where: 



tan Oi. = 



sin(fc) 
cos(fc) — g ' 



(82) 



(83) 



corresponding to the eigenenergies — ±ek with Sk — 
2\/l + 9 2 — 2gcos(k). It is easy to check that: 



sin(fc) 



2 1 + g 2 - 2< ? cos(fc)' 



(84) 



where we have used d g = ^§^dg h . 

For a linear quench protocol g(r) = 1 + A = 1 — vt, 
the imaginary time Schrodinger equation ([3| splits into 
a sum of independent differential equations: 



In the adiabatic limit, where only low momentum 
modes contribute to the excitations, we work with the 
linearized spectrum of the Ising model. Then the Hamil- 
tonian (81 1 simplifies to rlk = 2vt<j z +2ka x and Eqs. (85) 
and (1861) reduce to: 



cifc = —2vrbk — 2kbk, 
bk = —2ka,k + 2vrbk- 



(89) 
(90) 



dfc = 2[1 — vt — cos(fc)]afc — 2 sin(fc)6fc, (85) 
bk = — 2sin(k)a,k — 2[1 — vt — cos(k)]bk- (86) 

In the limit of t — > — oo we have 9 k — > and the eigen- 
states (82) simply become: 



(87) 



At the critical point g c = 1, corresponding to r — 0, we 
have 9k = - so that: 



i+>r° 



sin(^) 
cos(^) 



\r=0 
/fc 



cos(^) 

f — 
4 



These equations can be solved exactly analytically (we 
point out that this problem is the imaginary-time coun- 
terpart of the half Landau Zener (LZ) problem ana- 
lyzed in Refs. HU 04j |45|) . It is convenient to rescale 
the variables, r — > t/^/v, k — > gv 7 ^ an d differen- 
tiate both of these equations with respect to time. 
We then obtain a q — (4t 2 + Aq 2 - 2) a q = and b q — 
(4r 2 + 4<? 2 + 2) b q = 0. Each of these equations is of 
the type y — (Ax 2 + c)y = 0, which has two generic so- 
lutions: yi(x) = e~ x iFx (f + |, |, 2x 2 ) and y 2 (:r) = 

(-l)- 1 /4( 2a; ) e - 2 1 F 1 (f + |,|,2x 2 ), where i*i(a,M) 
is the confluent hypergeometric function. The generic 
solutions for Eqs. (p9| and (90) are therefore: 



o,(r) = cie- r ifi ( ^, -, 2r 2 ) + c*(-l)-V«(2r)e- r ^ 



(91) 



6g(r) = + \, i,2r 2 ) +c 4 (-l)- 1 /4 (2T)e -^ lFl (I + lf 3 j2t2 ^ 



(92) 



The coefficients c\ , c 2 , C3 , C4 are determined by the ini- 
tial conditions on the wave function at time To and 
two auxiliary conditions, e.g., a 9 | T= o = 6^(0) and 
6 9 | T=0 = — 2qa q (0) (continuity at t = 0). These last two 
equations give c 2 = (— l) 5 ^ 4 qc 3 , c 4 = (— l) 5 / 4 gci, while 
ci and C3 are set by the requirement that the system was 
in its ground state in the distant past: a g (r — > —00) = 1 
and 6 9 (tq — > —00) = 0. 



Using the expansion of the hypergeometric function 
when \z\ — > 00: 

TTi / 7 \ z a-b r (b) .-a T(6) n fl\ 

lFl (a,b,z) = ez f^ + i-z) f(b^aj + {z) 
we find that: 

q lV/2 + 1/2) 

C3 = ~7I r (g 2 /2 + i) Cl - (93) 
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Observable 


vL 2 < 1 


vL 2 > 1 


vL 2 > 1 




Exact sol. and APT 


Exact sol. 


APT 


M x 


[l/16]vL 2 


0.26</vL 


0.296^L 


Q 


[7<(3)/128tt> 2 L 3 


0.0265 vL 


0.0273 vL 


F 


[l/6144]u 2 L 4 


0.0276^^ 


0.0314^L 



TABLE I: Scaling of several observables in the ID trans- 
verse field Ising model after a linear quench in imaginary 
time [Eq. with quench velocity v and system size L. 

The second and third column show the asymptotics of the 
exact solutions for: the finite-size adiabatic limit, vL 2 <C 1, 
and thermodynamic adiabatic limit vL 2 2> 1. The last col- 
umn shows the scaling results within APT for vL 2 2> 1. For 
vL 2 <C 1 the perturbative (APT) and exact expressions are 
identical. 



The resulting probability of being in the excited state at 
the end of the evolution for r = is found by overlapping 
the final wave function with the excited state at t = 0: 



+ 



T=0 



1 

71 



(94) 



1. Excess energy Q 

We consider the total excess energy of the system (en- 
ergy above the instantaneous ground state energy): 



(H) - (H) Q . 



(97) 



The scaling predictions in Eq. ( p5| ) apply to this case 
when the observable M 7 = (i/j\J\A^\i/;) is associated with 



the Hamiltonian operator: M.„ 



H, according to our 
(22 1. In this case, 



definition of generalized force in Eq 
since we are dealing with a diagonal operator, the geo- 
metric tensor is identically zero [see Eqs. (26) and ( 34 )] 



and the response is quadratic and, in this case of a quench 
of a single parameter, is proportional to a single compo- 



nentn^ja (Eq.[30j 



y> \(0\dxH\n) 



(98) 



The response to the quench, both in real and imaginary 
time, is: 



Therefore we have: 



-PCX 



(<z) = 



1 |q g (0) + 6 g (0)| ; 

2M0)p 
i 



HM0)l 2 

q r( g 2 /2+i/2; 



V2 r(g2/2+l) 



1 



q r(g 2 /2+l/2) 
v/2 r(«a/2+l) 



2 - 



(95) 



where we point out that, since i-Fi(a, 6,0) = 1 for any 
a and b, we have a 9 (0) = c\ and b q (0) — C3. Restoring 
the dependence on v using the substitution q — > k/y/v 
we obtain: 



p£f(M) 



r c+i — 



7 + J 



(96) 



(99) 



From the scaling dimension of the susceptibility n^^^ it 
is possible to extract the scaling behavior of Q, as already 
done in Ref. [3], and to take into account the finiteness of 
the system, as discussed in Section \nj\ without knowing 
the details of the model. The expected scaling behavior is 
stated in Table [TJ By applying APT to the specific Ising 
Hamiltonian under investigation, we could also extract 
the numerical prefactors of the scaling, as we will explain 
in the next section. 

From the exact solution presented in section |V A[ we 
evaluate the total excess energy Q at the final critical 
point in the scaling limit as: 



Q 



k>0 



eIp l J{Kv) 



k>0 



4* p£f (*,«). 



(100) 



Furthermore we note that to correctly define the fi- 
nal amplitudes a+._(r = 0) on the |+) and |— ) eigen- 
state, we need to properly normalized the coefficients in 



Eqs. (91) and (92) to satisfy the condition in Eq. (15). 



In the limit vL 2 3> 1, we convert the sum into an integral 
to find 



Q = / dqqp^{q) = 0.0265L?;. 

n Jo 



(101) 



B. Observables 

In the following we present some observables that 
describe the response of the system to the quench in 
Eq. ( 79 ) . Their scaling behavior is derived from the ex- 
act solution using the excitation probability in Eq. ( 96 1 



The results are summarized in Table |TJ where they are 
also compared with the correspondent scaling found by 
adiabatic perturbation theory (APT). 



In the limit vL <C 1: p^f(q,v) ~ and the total 

excess energy becomes: 



4^2 
"64" 



OO 1 

V 1 

^ \U2m + V) 



, 2r3 7C(3) 



,0 L i 



f (2m + l)] 3 V L 128tt 3 ' 



(102) 



where we used anti-periodic boundary conditions for 
fermions which map to periodic boundary conditions for 
spins (see Ref [46 for details). 
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2. Log-fidelity 

The logarithm of the fidelity F = - ln(| (V>(0)|0) | 2 ), 
in the perturbative regime that we are considering, ac- 
cording to Eq. (231, can be approximated as: F ps 
J2 n ^o \ a n\ 2 - Therefore the scaling of the log-fidelity can 
be extracted from that of the generalized force corre- 
sponding to the identity operator: JvL, = —I, as it is 



easy to see from our definition in Eq. ( 22 ) . Based on the 



same reasoning as above for the excess energy, the linear 
response for this observable vanishes, and the coefficient 
for the quadratic scaling is: 



FXX 



£ 

n^O 



\(0\dxH\n)\* 



(103) 



From the exact solution (see also Ref. [47] ) we can cal- 
culate F using: 



F = -£>(l-j& z (M)). 

fc>0 



(104) 



Then in the limit vL 2 ^> 1, transforming the sum to 
an integral we immediately find F = 0.0276\/vL, while 



in the opposite limit F 



64 



fil 14 



v 2 L\ 



(2m+l) 4 



3. Transverse Magnetization 

Finally we study the transverse (excess) magnetization 
M S = 5>?>-£< KI >- (105) 



According to the definition in Eq. (22 1, M x corresponds 



to the expectation value of the observable M. x = —d\H, 
i.e., the generalized force with respect to the coupling 
constant that is quenched in time, which in our case is 
g(r) = 1 + A(t). Therefore, from Eq. (25) we expect the 
scaling: 



M x fa -2v x g\ t \, 



(106) 



where it should be noted that v\ = d T X\ = —v. We 
extract the value of M x from the exact solution a follows: 
we evaluate for each momentum k the expectation value 
at the end of the process (r = 0): (a z (k)) T=0 = — (cf — 
c i)/( c i + c 1)j an d the sum over all momenta: 



M x 



r( 



L ) 2v L \2v 1> 



r( 



i) 



■-r(— + 

2v V 2v ~ 2 



l\2 



(107) 



Evaluating this expression in the limit vL? 3> 1 we find 
M x = 0.264-^27, while in the opposite limit M x pa 



4\irJ " ^im=0 (2m+l) 2 

These results from the exact solution are compared 
m Table [T] with the ones from APT. The agreement is 



very good, we will comment on this in more detail in the 
following section. 

We point out that if instead we would perform a quench 
changing J(r) = 1 — A = 1 + vt, with g = 1 (as we did 
with QMG in Ref.|3]), the correspondent generalized force 
Ai\ is now found with A = —J. It corresponds to the 
observable: 



E, 



-J 



(ij) 



£< 

(ij) 



(108) 



which is the excess interaction energy, or z-energy. 
Therefore, the observables E z and M x have the same 
scaling behavior respectively for a quench of the coupling 
J and of the coupling g. When we deal with the QMC 
simulations and the numerical solution of the differential 
equations, it is in practice more convenient to perform 
a quench of the J coupling, therefore in the following, 
when presenting the data we will use the observable E z . 



C. Adiabatic perturbation theory for the 
transverse field Ising model 

The exact solution for the quench dynamics in imagi- 
nary time of the ID transverse-field Ising model provides 
a good opportunity to test the APT method presented in 
Sec. [n] The APT analysis in this case is very accurate, 
agreeing very well with the exact solution, as shown in 
Table [TJ The bas ic ingredient of the APT is the transi- 
tion amplitude a n (rf = 0), which in terms of the tuning 
parameter A = —vt becomes [3]: 



a n (0) w / d\(n\d x \0)e X p 



dX' 



(E n (X') - 8 m (X') 



(109) 

For the transverse-field Ising model the lowest excitations 
correspond to flipping the effective spin from \—)k to \+)k 
(corresponding to exciting two Bogoliubov's fermions 
with opposite momenta) characterized by the matrix ele- 
ment in Eq. ( |84[ | and the energy difference 2eu- Therefore 
the transition amplitude from the ground to the excited 
state is: 



afc(0) 



dX 



sin(fc) 



A 2 + 2(A+l)(l-cos(fc)) 



X exp"* So dAyA 2 +2(A+l)(l-cos(fc)) _ 



This expression simplifies in the slow limit, where we can 
use the linearized spectrum: 



afc(0) ^/ dA A^ eXP 



dXWk 2 + A' ; 



(111) 
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Using these transition amplitudes instead of the exact 
expressions found in the previous section we recover the 
last column in Table [J 

In Fig. [3] we plot the excess energy Q for different sys- 
tem sizes as a function of vL 2 computed from the exact 
solutions of the Schrodinger equation [solving numeri- 
cally Eqs. (89) and (90 1]. In the same plot we also show 



the APT results for the infinite size. The agreement be- 
tween the two methods is excellent when sufficiently large 
systems are used for the exact solution. 

In most situations the exact solution for a time- 
dependent problem is not readily available. The good 
agreement we found here suggests that one can make 
many qualitative and even quantitative statements about 
the dynamics using the adiabatic perturbation theory, 
which only requires the integration of static quantities. 



10 



10 l 



10 



QL 




10 



FIG. 3: (Color online) Results from the exact solution of the 
Schrodinger equation. The excess energy Q following an imag- 
inary time quench for different system sizes L = 8, 16, 32, 64 
is shown as a function of vL 2 , with v varying between 128 and 
10 -4 ). The data collapse is evident, the splitting of the curves 
for large v L 2 is due to finite size effects. The points overlap 
well with the predictions from APT based on the linearized 
spectrum, shown with dashed and dot-dashed straight lines 
for the high- and low- velocity regimes. The fitted lines have 
slope 2 for vL 2 <C 1 and 1 for vL 2 2> 1, in agreement with 
the APT scaling in Table [I] 



VI. OBSERVABLES IN REAL AND 
IMAGINARY TIME QUENCHES OF THE ISING 
MODEL 

According to the results presented in the previous sec- 
tions, it follows that the imaginary-time dynamical pro- 
tocols give very similar results as the real-time protocols 
considered earlier |14) . This implies, in particular, that 
imaginary-time quantum evolution with a time evolv- 
ing Hamiltonian can be used for simulations of real-time 
non-equilibrium dynamics, including, e.g., realizations of 



the quantum Kibble-Zurek mechanism [2S]- However, 
there are still also important differences between the two 
types of dynamics. Firstly, the imaginary-time evolu- 
tion clearly breaks time-reversal symmetry and this in- 
troduces a strong asymmetry between the initial and fi- 
nal times of the evolution. Thus, in real-time dynamics, 
the system is constantly excited during the evolution and 
these excitations propagate in time. On the contrary, 
during the imaginary-time evolution the system always 
relaxes toward the ground state, and the effects of non- 
adiabaticity are visible only when approaching the final 
state, i.e., the critical point. For example, in real-time 
evolution the ground state fidelity and the diagonal en- 
tropy [48] (which are observable independent measures of 
non-adiabaticity) are identical for the time-reversed pro- 
tocols. In particular, the degree of non-adiabaticity is the 
same if one considers protocols which start or end at the 
quantum critical point. In imaginary-time evolution this 
is not the case. If one passes a singularity, like a critical 
point, in a real-time process, then it will always result in 
non-analyticities in various observables (the defect den- 
sity in the case of Kibble-Zurek mechanism is an example 
of this). In imaginary-time evolution the singularities in 
the observables will show up only if one ends the process 
at this singularity or in its close vicinity. 

In this section we analyze closely the behavior of the 
observables in the case of quenches of the transverse-field 
Ising model, comparing the exact solutions for the real- 
and imaginary-time cases. For simplicity we consider 
here the protocol already analyzed in Ref. [3J, where we 
fix g = 1 and ramp J linearly in time to end at the critical 
point, i.e., J = 1 + vt and J = 1 + vt, in imaginary and 
real time, respectively, with the final time: r/ = tf — 0. 
Then the scalings of the excess energy, fidelity and mag- 
netization in real and imaginary times look nearly identi- 
cal. Since real-time evolution in this model was analyzed 
earlier in different papers [TJJfSS] we will omit the details 
of the calculation and only present the final results. We 
note that in Ref. [14] we analyzed a linear quench where 
one starts at the quantum critical point. Because of the 
symmetry of the transition probabilities with respect to 
time reversal the analysis applies as well to the process 
we are interested in here, where one ends the quench at 
the quantum-critical point. The only subtlety appears 
in the analysis of the x-magnetization, which is an off- 
diagonal observable and which depends on the phase of 
the transition amplitude. We will comment on this sub- 
tlety below. 

In Table [TT] we present the comparison of the scaling 
of several observables for linear quenches to the QCP in 
real and imaginary time, obtained from the exact solu- 
tion of the transverse-field Ising model (see Sec. 
for the imaginary-time case and Refs. [131 [35] for 



VA 



the 



real-time case) . In Figures |4j [5j and [6] we plot the corre- 
sponding quantities obtained by solving numerically the 
Schrodinger equation. The definition of the observables 
was given in the previous Section [V] As mentioned be- 
fore, since we are quenching J (and not g, to be consistent 
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TABLE II: Results from the exact solution of the ID 
transverse-field Ising model: scaling forms for the magneti- 
zation, the excess energy, and log-fidelity with the quench 
velocity v and the system size L in real and imaginary time 
for different regimes. 



with the protocol used in QMC simulation presented in 
following sections) , the observable which gives the fidelity 
susceptibility in the linear response is the excess interac- 
tion energy E z [see Eq. (108)]. It is expected to scale in 
the same way as M x . 
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FIG. 4: (Color online) Comparison of real- and imaginary- 
time dynamic scaling of the excess energy QL (left) and the 
log-fidelity F (right) for system size L = 16 (top) and L = 64 
(bottom). 



Overall we find very good agreement and almost iden- 
tical behavior between the imaginary and the real time 
cases. A more careful analysis is nevertheless necessary. 
For the diagonal observables, the excess energy Q and 
the log-fidelity F, the scaling behaviors are the same in 
real and imaginary time and in agreement with the APT 
predictions presented in the previous section. In partic- 
ular, in the limit vL 2 <C 1 even the prefactors coincide — 
indeed, in this limit the analytic expression are identical. 
In the opposite regime vL 2 ^> 1 the prefactors are slightly 
different. In this limit, the real-time dynamics presents 
a more oscillating behavior: see for instance the plots of 
the excess energy in Fig. [5] A similar behavior was also 
observed in Refs. QUI El]. 

The case of the excess x-energy or magnetization along 
the x-direction [as defined in Eq. ( 105 )] requires more 
attention. Indeed this quantity, as mentioned before, 
corresponds to the generalized force with respect to the 
coupling A that drives the dynamics. Working out the 
asymptotic scaling behavior from the scaling dimension 
in the limit of vL 2 » 1 we find M x ~ y/vL in both 
real and imaginary times, according to Eqs. (24 1 and 



(25). Concerning the limit vL <C 1, in imaginary time, 



from the exact solution [see Eq. ( 107 )] we know that 
M x ~ vL 2 . In the real-time case, from analyzing the 
exact solution we can infer that the behavior for small 

^3 

vL 2 is non analytic, decaying exponentially as ~ e~ . 
Such behavior is visible in the plot in Fig. [6j for large 
values of vL 2 (but not too large, as finite-size effects also 
are apparent) the slopes of the real- and imaginary-time 
functions are the same, the data being shifted by a fac- 
tor of 2 according to the predictions. For vL 2 <C 1 the 
imaginary-time function decays analytically with slope 1 
as expected, while in the real-time case there is a more 
rapid drop reflecting the non-analyticity of the function. 
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FIG. 5: (Color online) Data-collapse plot for the excess en- 
ergy Q (top) and the log fidelity F (bottom) based on real- 
(left) and imaginary-time (right) dynamics for different sys- 
tem sizes. In the regime of large vL 2 the splitting of the 
curves is due to finite-size effects. The real-time case show 
more oscillating behavior than the imaginary-time case. 



VII. APPLICATION: DETECTION OF 
QUANTUM-CRITICAL POINTS 

A useful application of the universal scaling presented 
in the previous sections can be found by considering 
quenches in either real or imaginary time that sweep 
through the QCP and end at different final amplitudes 
Xf =/= X c . In this discussed in Sec. IIIIl the length 

scale £a(A) ~ |A— A c | _!/ is not diverging anymore and par- 
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FIG. 6: (Color online) The excess interaction energy E z 
(which has the same scaling as the magnetization M x ) in real- 
and imaginary-time dynamics for system sizes L — 32 (black) 
and L — 64 (red). For vL 2 3> 1 the slopes are the same, 
except for a shift due to a different pref actor. The splitting 
of the curves for large vL 2 is due to finite-size effects. For 
vL 2 <C 1 the real-time case changes drastically, decaying to 
zero exponentially. 



ticipates in the scaling behavior along with ~ \v\ *"+! 
and the system size L. For a generalized force M 7 we then 
expect a scaling behavior as the one given in Eq. (66): 



M 7 (X,v) = coBSb + L-i*UL/Z x (\),L/Z v ) 



const 



+ Lv hy v i/(zv+i)> VL J' 



(112) 



as was already suggested in Ref. If we perform 

several quenches changing the final amplitude A/ and 

plot for each of them the rescaled quantity M 1 /v 1 + uz , 
we expect all the curves (asymptotically for sufficiently 
large L) to cross at the location of the QCP, since when 
Xf = A c the rescaled quantity does not depend on v 
anymore. We have performed such an analysis for an 
imaginary-time quench of the form J(r) = 1 — A = 1+vt, 
and the correspondent real-time one (replacing r with t) . 
Sweeping across the critical point, i.e., starting from a 
negative A and ending at some positive value, and look- 
ing at the z-energy E z as defined above in Eq. (108), we 
expect: 



E m 



,vL' 



(113) 



In Fig. [7] we show the results for imaginary- (top graph) 
and real-time (bottom graph) quenches based on the nu- 
merical solution of the Schrodinger equation. The ex- 
pected behavior is confirmed in both cases, all the lines 
cross around J = 1 where the QCP is located. There- 
fore, through this type of analysis it would be possible to 
locate the position in the parameter space of a QCP of 
a system that cannot be solved exactly and of which the 
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FIG. 7: (Color online) Real-time (bottom graph) and 
imaginary-time (top panel) results for quenches ending at 
variable finite amplitude Jf. The rescaled z-energy E z /y/v 
is shown for different quench velocity v and system size L, 
with the product vL 2 fixed at 512. In both cases the curves 
cross around the location of the quantum critical point J = 1. 
The inset shows the same data with the :r-axis rescaled ap- 
propriately to achieve the data collapse. 



position of the critical point is not known. The univer- 
sal behavior described by Eq. (113) is furthermore con- 



firmed by the collapse of the data when plotting versus 
the rescaled quantity (1 — Jf)y/v: see the insets of Fig. [7] 
As expected, a similar collapse is also observed if we per- 
form a quench that does not reach the QCP but ends just 
before it, as we show by the data in Fig. [8j These results 
were obtained by numerical simulation of the ID Ising 
model with the NEQMC method introduced in Ref. [3] 
(and discussed also above in Sec. IV), using it to perform 
a linear quench ending at different values of the final am- 
plitude Jf and approaching J c = 1. 
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FIG. 8: QMC results: data collapse for the z-energy E z for a 
linear quench of the coupling J ending at different values J/ 
just before the QCP. 



FIG. 9: Universal relaxation dynamics in imaginary time after 
a linear quench. 



VIII. UNIVERSAL RELAXATION TO 
EQUILIBRIUM AFTER A QUENCH 

Up to this point we have been concerned with the scal- 
ing of the observables right at the end of a quench, in 
particular at the final time, when the Hamiltonian of the 
system has reached the QCP \f = A c or, as in the pre- 
vious section, some other amplitude A/ =/= A c . Here we 
want to address the scaling behavior that follows sub- 
sequent to the quench, when the system starts relaxing 
governed by a fixed Hamiltonian. Therefore we let the 
system evolve after the end of a quench with the fixed 
final Hamiltonian for a variable length of time, that we 
call t R , and we look at the behavior as a function of t R , 
that wc call the relaxation time. Based on the scaling 
arguments that have lead us to Eq. (66), we argue that, 



if we let the system evolve after the quench for a time 
t R , then we expect [for a generic r-th power quench as in 
Eq. |53b]: 



M 7 (A,u) ~ (t R v^,L z /t R , \X f - A c | 



" v t R ) . 
(H4) 

This means that the relaxation time itself comes into 
play in the universal scaling behavior as an additional 
"length" scale to be compared to the other characteristic 
lengths of the system. For sudden quenches this conjec- 
ture was recently suggested and tested in Ref. [35]. If 
the system size is large enough, for instance, we expect 
the quantity t R v z^+i to be the rescaled variable that 
characterize the universal relaxation after quenches with 
different velocities. As before, we have checked this be- 
havior for the z-componcnt excess energy E z from the 



exact solutions [solving numerical Eqs. (89 1 and (90)] for 
a linear quench, see Fig. [9] and with NEQMC for a sud- 
den quench see Fig. [lO] 
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FIG. 10: QMC data: relaxation dynamics of the interaction 
energy E z in imaginary time after a sudden quench. 



IX. SUMMARY AND CONCLUSIONS 

We have presented an overview of general aspects of 
non-adiabatic response of physical observables to slowly 
changing parameters, both in imaginary (Euclidean) and 
real time. There are many similarities between the 
imaginary- and real-time response, which we demon- 
strated by calculating the leading first- and second- 
order non-adiabatic corrections of physical observables. 
We identified the corresponding susceptibilities and ex- 
pressed them through the non-equal time correlation 
functions. In particular, we extended the traditional 
Kubo response theory to describe the response of systems 
to perturbations which are slow but can be arbitrarily 
large in amplitude. The components of the geometric 
tensor (the metric tensor and the Berry curvature) natu- 
rally emerge as response functions of physical observables 
to the quench velocity. 
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Both real- and imaginary-time dynamics near contin- 
uous phase transitions can be used to analyze univer- 
sal non-adiabatic response of the system and extract 
static and dynamic critical exponents, using a general- 
ized non-equilibrium scaling theory, which we also fur- 
ther elaborated here. Importantly, imaginary-time dy- 
namics is amenable to powerful Monte Carlo simulation 
methods. We briefly reviewed two different QMC algo- 
rithms which directly implement quantum dynamics for 
interacting systems. They have the same range of prac- 
tical applicability (avoidability of sign problems) as con- 
ventional equilibrium finite-temperature or ground-state 
projection methods. 

We illustrated the utility of the general theoretical for- 
malism using the particular example of the transverse- 
field Ising model in one dimension. Using both exact 
treatments (through the standard mapping to fermions) 
and QMC simulations, we found that imaginary- and 
real-time dynamical responses indeed are very similar 
near the critical point, for all physical observables ex- 
amined. We also found excellent agreement with pre- 
dictions of adiabatic perturbation theory. We illustrated 
how one can use the non-equilibrium finite size scaling 
to accurately extract the transition point and the critical 
exponents. 



The ideas presented in this article have many potential 
applications, including (i) analysing universal dynami- 
cal response near quantum phase transitions with un- 
known dynamical exponent, e.g., in disordered systems; 
(ii) applying QMC methods to implement imaginary- 
time quantum annealing and (using the similarity of non- 
adiabatic response in real and imaginary times) making 
predictions concerning real-time quantum annealing pro- 
tocols; (iii) using non-adiabatic response of physical ob- 
servables to directly extract the Berry curvature and the 
metric tensor (including the fidelity susceptibility) either 
experimentally or numerically. 
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